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Experimental realizations of disorder in optical lattices generate a distribution of the Bose- 
Hubbard (BH) parameters, like on-site potentials, hopping strengths, and interaction energies. We 
analyze this distribution for bosons in a 2d bichromatic quasi-periodic potential by determining 
the generalized Wannier functions and calculating the corresponding BH parameters. Using a local 
mean-field cluster analysis, we study the effect of the corresponding disorder on the phase diagrams. 

We find a substantial amount of disorder in the hopping strengths, which produces strong deviations 
from the phase diagram of the disordered BH model with purely random on-site potentials. 


PACS numbers: 

I. INTRODUCTION 


Bosons in a regular optical lattice, which can be de¬ 
scribed by the BH model, show a transition from an in¬ 
compressible, non-coherent Mott insulator (MI) phase to 
a compressible, coherent superfluid (SF) phase, which 
has been shown experimentally [1] as well as theoreti¬ 
cally [2, 3]. In the disordered case the Bose-glass (BG) 
phase occurs between the MI and the SF, which is com¬ 
pressible, but not coherent and whose existence has been 
shown experimentally [4, 5]. 

Theoretically, disordered BH systems are modeled with 
isolated on-site box disorder, while the other parameters 
have precise values. Numerical methods like quantum 
Monte Carlo (QMC) methods [6-14] and density-matrix 
renormalization group techniques [15-19] have been ap¬ 
plied to this system in order to study the phase diagram. 
Based on the mean-field approximation [20] various nu¬ 
merical techniques, such as local mean-field (LMF) the¬ 
ory [21-23], stochastic mean-field theory [24, 25], and 
LMF cluster analysis [26] have been proposed. The phase 
diagram of the 2d disordered BH model with random on¬ 
site potentials predicted by the LMF cluster analysis [26], 
which interprets the BG-SF transition as a percolation of 
SF regions, agrees very well with the one obtained from 
QMC methods [9]. The authors of [27] apply the perco¬ 
lation criterion for their studies on the BH model with 
random impurities. 

Disorder in other BH parameters, like the tunneling rates 
or the inter-particle interactions, has not been studied as 
intensively as disorder in the on-site potentials. Disorder 
in the inter-particle interaction can be realized experi¬ 
mentally near a Feshbach resonance [28, 29]. A uniform 
distribution of this parameter has been studied in [28, 30]. 
Disorder of the tunneling rates has been studied within 
SMF theory [25] for a uniform distribution. Other works 
on disordered tunneling rates focus on bimodal distribu¬ 
tions. 

Experimentally, disorder can be introduced either by a 
speckle field [31] or by a bichromatic quasi-periodic lat¬ 
tice produced by two lasers with incommensurate wave 
lengths [4, 32]. While the distributions of the BH pa¬ 


rameters resulting from speckle fields have been studied 
in [31, 33], in this paper we focus on bosons in two- 
dimensional bichromatic quasi-periodic potentials and 
analyze the distribution of all of the BH parameters, i.e. 
on-site potential, hopping strength and interaction ener¬ 
gies, by determining the generalized Wannier functions 
and calculating the corresponding BH parameters. 

This paper is organized as follows: First, we discuss the 
BH model and the LMF cluster analysis in section II. In 
section III we determine the phase diagrams of the disor¬ 
dered BH model for three different cases of uncorrelated 
disorder for comparison: exclusively random on-site po¬ 
tentials, exclusively random hopping strengths and exclu¬ 
sively random interaction energies. In section IV the BH 
parameters are calculated for a two-dimensional bichro¬ 
matic quasi-periodic potential and their distributions are 
characterized. These are finally used in section V to de¬ 
termine the phase diagrams. The paper concludes with 
a discussion. 


II. THE BOSE-HUBBARD MODEL 

The BH Hamiltonian describing Bosons in an optical 
lattice [34] is given by 

# = + j TAG -1 )-(!) 

i i < i,j) 

The operator h i = a\a i is the particle number operator 
of bosons on site i, which are annihilated and created by 
the operators a i and a\. The site index i = 1,..., M, 
where M = L 2 is the number of sites in a L x L 2d lat¬ 
tice, represents a tuple of spatial coordinates (i x ,i z ) with 
i X:Z = 1,..., L. The chemical potential is denoted by /q 
the inter particle repulsion by V and the tunneling rate 
by J. The last sum runs over all four (Z — 4) nearest 
neighbor pairs (q j) of the lattice. 

In the ordered case, the fi/U-JZ/U phase diagram dis¬ 
plays the well-known Mott-lobes with boundaries given 
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in LMF theory by [2, 3] 

IM±{J,Z,U,n) = - l -(JZ-U(2n-l)) 



Inside the Mott-lobes the particle number is fixed to an 
integer value n. They are aligned along the fi- axis and the 
particle number increases from lobe to lobe with grow¬ 
ing [i. In this region particle tunneling is prohibited due 
to the existence of an energy gap in the particle excita¬ 
tion spectrum. The Mott-lobes are surrounded by the 
SF phase, where particle tunneling is favorable and co¬ 
herence grows with increasing tunneling rate. 

In the case of on-site disorder, the BG phase occurs be¬ 
tween the MI and SF phase. In this phase the system is 
compressible, but not SF. In the present paper we will 
not restrict to on-site disorder, but rather study the influ¬ 
ence of disorder on all BH parameters. Therefore, all BH 
parameters, except the chemical potential /q which is a 
global parameter fixing the particle number, now become 
site dependent: 

H = ^2{ei- n) fit + y 5^ («i - i) - 53 JiA^y 

i i (i,j) 

(3) 

The SF parameter 

A = (di) (4) 

is defined as the expectation value of the annihilation 
operator and can be chosen to be real because of the 
U (1) "symmetry of the BH-Hamiltonian. With the help 
of the LMF approximation [20] 

a A « AA + djipi - Ai’j, (5) 

the Hamiltonian can be transformed into a sum of local 
Hamiltonians H = 

H, = (e» - m) A + y hi (y Al) ~ V i (a» + 4 ~ A) , 

(6) 

which are effectively coupled by a local hopping rate 
r]i := JT JijAij'ipj , which depends on the local SF pa¬ 
rameter of the neighboring sites with Aij = 1 for % and j 
being nearest neighbors on the square lattice with peri¬ 
odic boundary conditions and zero otherwise. One should 
note that mean field approximations as the one used 
here neglect spatial correlations of quantum fluctuations, 
which is less severe in high space dimensions. Actually it 
is exact in infinite dimensions and reproduces the exact 
critical behavior already above the upper critical dimen¬ 
sion, which is unknown for the system we consider here 
but certainly larger than 4. Consequently, we expect the 
approximation to be critically inaccurate for Id systems 
and therefore, focus on two space dimensions. We would 
expect the predictions of our LMF analysis to be even 


more accurate for three dimensions, but 3d systems are 
computationally much more demanding. 

Since the Hamiltonian itself depends on the SF parameter 
[22], the self-consistency equation (4) is solved recursively 
in order to find the ground state of the LMF Hamiltonian 
(6). In turn all local parameters of interest, especially 
the local particle number rq = (ajcq), can be computed. 
In order to determine the phase boundaries, several steps 
are needed [26]: First, we define and identify so-called MI 
and SF sites. Sites with an integer number of particles 
are called MI sites. Note that the LMF approximation 
neglects quantum fluctuations and predicts an integer ex¬ 
pectation value for the local particle number for some site 
even in the presence of disorder in the BH parameters. 
Sites with non-integer expectation value of the local par¬ 
ticle number are denoted as SF sites. In a second step 
all three occurring phase are identified: In the MI phase 
the system only consists of MI sites and no SF site oc¬ 
curs. The BG is characterized by a mixture of MI and 
SF sites, more precisely by isolated clusters (islands) of 
SF sites within a sea of MI sites. In [26] we discussed 
how this definition is physically plausible regarding the 
conventional hallmarks of the BG phase, namely the lack 
of coherence and a gapless spectrum. Since the particle 
number fluctuates within the isolated SF clusters only 
the sites within a single cluster can be phase coherent, 
sites in different clusters are not, which establishes the 
lack of macroscopic phase coherence. Moreover, the BG 
phase is actually the Griffiths phase of the BH model 
(see [26] for a discussion), where the isolated SF clusters 
behave like finite systems within the SF phase and there¬ 
fore, have a very small gap that decreases quickly with 
the size of the SF cluster. Since SF clusters can be arbi¬ 
trarily large there is no lower bound for the gap and thus 
the BG phase is gapless. While approaching the BG-SF 
transition the regions with SF site grow, form connected 
clusters, which finally percolate. The percolation of the 
SF sites marks the transition to the SF phase. We call 
this scheme the LMF cluster analysis approach, which 
is described and discussed in [26]. In contrast to other 
LMF approaches we show, that the LMF cluster analy¬ 
sis reproduces the phase diagram predicted by quantum 
Monte Carlo [6] with excellent accuracy. 


III. BOX DISTRIBUTED DISORDER AND 
DIFFERENT SCENARIOS 

Before we focus on the bichromatic quasi-periodic 
potential, we analyze the effect of uncorrelated disor¬ 
der in each of the BH parameter q, Jq- and Ui sep¬ 
arately. Here we determine, with the LMF cluster 
analysis, the phase diagram for a uniform distribu¬ 
tion p{a) = 0 (A a /2 — |a|) /A a for each BH parameter 
a = ei, Jij , Ui separately. While all three phases can be 
found in each disorder scenario, we find substantial dif¬ 
ferences in the phase diagrams. 
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A. Disorder in e 

The most common disorder scenario is diagonal dis¬ 
order introduced by site-dependent local on-site en¬ 
ergies q, which are drawn from a box distribution 
p(ei) = 0 (A e /2 — 16^|) /A e . This has been widely stud¬ 
ied via quantum Monte-Carlo methods [6-8, 12], mean- 
field techniques [22, 24, 26] and analytic approaches 
[9,35]. 

Figure 1 shows the phase diagram for different disor- 



FIG. 1: Phase diagram for box distributed disordered on-site 
energies ei with the disorder strength A e /U — 0.35, 0.6, 1.5. 

der strengths resulting from the LMF-Cluster-Analysis 
[26]. The Mott-lobes with a fixed particle number 
of n = {n/) av extend from /i_ = (n — 1) V + A e /2 to 
fi + = nU — A e /2, see [2]. Thus, all Mott-lobes simul¬ 
taneously disappear at a critical disorder strength of 
A C JU = 1. The Mott-lobes are surrounded by a BG 
region. For larger tunneling rates a phase transition to 
the SF regime occurs. While all three phases appear for 
A C e /U < 1, in the strong disorder limit only the BG in the 
small tunneling and the SF phase in the high tunneling 
regime survive. 

B. Disorder in J 

The influence of disordered tunneling rates has mainly 
been studied for bimodal distributions, where two values 
of the tunneling rate are chosen and distributed randomly 
among the lattice [10, 25, 36, 37] leading to fundamen¬ 
tally different phase diagrams compared to the one dis¬ 
cussed here. In contrast, here we focus on a general ap¬ 
proach, where the local tunneling rates are uniformly dis¬ 
tributed according to p(Jij) = J + 0 (Aj/2 — \Jij\) / A j 
symmetric around a fixed value J. 

In Figure 2 the phase diagrams resulting from LMF- 
Cluster-Analysis in dependence of the chemical potential 
ji/U and the mean tunneling rate JZ/U are shown. It is 
important to notice that here the width of the disorder 



FIG. 2: Phase diagram for box distributed disordered 
tunneling rates Ji with the disorder strength A j/U = 
0.035, 0.06, 0.095. The black line is the perturbative result 
for MI-SF transition in the ordered system given by (2). 


distribution is one order of magnitude smaller than for 
on-site disorder. Here two new unique features occur in 
the phase diagram: Firstly, the BG regions are separated 
into individual regions by SF regions reaching down to 
JZ/U -A 0. Secondly, the distance between the Mott- 
lobes increases with their number n and the number of 
Mott-lobes is finite. For an intuitive explanation, we re¬ 
call the perturbative result of the ordered case, given in 
equation (2), which is the limiting result for vanishing 
disorder strength A j -A 0. For a small but non-zero 
disorder strength A j, the first SF regions in the lattice 
occur at sites with J + Aj/2. Consequently, the bound¬ 
ary of the MI region is shifted to the left by — Aj/2 with 
respect to ordered case with tunneling rate J. Since in 
the ordered case the tip of the Mott lobes (at tunnel¬ 
ing rate J m ax(^)) decreases with increasing // as l//i the 
shifted Mott lobes for tunneling disorder strength A j 
must disappear when 1/p gets smaller than A j resulting 
in only a finite number of Mott lobes. With other words: 
Even at a vanishing average tunneling rate (J = 0) a 
non-vanishing disorder in the tunneling rate produces an 
increasing fraction of SF sites when /i is sufficiently large. 
To put it quantitatively: Along the /r/t/-axis (J = 
0) the MI-lobes exist between /x_ (Aj/2, Z, ?7, n) and 
/i + (Aj/2, Z, U, n), with fi± given by equation (2). 
For fixed disorder strength A j the height of the 

Mott-lobes is given by A^ 1 = 2 yj^ — l) 2 — 

which decreases with n. The width becomes zero for 
n c Ll = \ (%cT — l) a7z’ which means that only a finite 
number n/f 1 of Mott lobes exist. As a consequence, the 
Mott-lobes disappear one after the other for increasing 
disorder strength A j. The last Mott-lobe (n = 1) disap¬ 
pears at ^ 0.0858. This is different from the 

on-site disorder case, where all vanish at the same critical 
disorder strength. In Figure 2 a) (A j/U = 0.035) three 
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Mott-lobes exist, two of which are visible, while in b) 
(A j/U = 0.06) only one remains. In the last diagram 
(A j/U = 0.095) no Mott-lobe exists, as the critical dis¬ 
order strength is exceeded. 

The Mott-lobes are surrounded by the BG phase. As a 
new feature in comparison to on-site disorder case we find 
disconnected BG regions between fi_ (Aj/2, Z = 1, t/, n) 
and /i + (Aj/2,Z = 1,17,n), with fi± given by equation 
(2), which are separated from each other by the SF region 
in the vicinity of integer values of fi/U. The fact that the 
SF region survives in the limit J —>• 0, is a unique feature 
of tunneling disorder. The width of the BG regions along 

the n/U-axis is given by A^ G = 2 \J\ (f % - l) 2 - f %n, 
which also decreases for growing n. The number of 
BG regions is given by n^ G = \ — l) Even 

though the BG regions survive for even higher disorder 
strength than the MI-lobes, they analogously disappear 
one after the other and finally disappear completely at 
= 2 (3 — 2\/2) « 0.3431. The SF phase exists for in¬ 
finitesimal small tunneling rates between these BG re¬ 
gions. At the ends of the BG regions narrowing tips 
occur, which are located along the line of mean integer 
filling but finally end in the SF region. 


C. Disorder in U 

Disorder in the inter-particle interaction can be real¬ 
ized near the Feshbach resonance [28, 29] and a uniform 
distribution of this parameter has been studied in 
[28, 30]. The phase diagrams for this case resulting from 
LMF cluster analysis are shown in Figure 3 for increasing 
disorder strength, where U is the mean value of the 
disorder distribution p(Ui) = U + 0 (A[//2 — \Ui\) /A jj. 
Analogously to the disordered tunneling case we find 
a finite number of Mott-lobes. Intuitively this can be 
understood by recalling the MI boundaries (2) of the 
ordered case, as we have already discussed for tunneling 
disorder. For small tunneling rates, the first SF sites 
occur, where the tunneling rate J overcomes the reduced 
inter-particle interaction U — A^y/2. Thus, the Mott- 
lobes shrink all by the same amount for fixed disorder 
strength Au. Therefore, the smallest Mott-lobes of the 
ordered system disappear leading to a finite number 
of Mott-lobes. Along the n/U— axis they extend from 
fi- = (n — 1) (U + to /i + = n(U — 4y) and they 
disappear at a critical disorder strength A' u/U = 
where and meet [28]. As the critical disorder 
strength Afj depends on the number n of the specific 
Mott-lobe, they vanish one after another, until for 
A \jjU = 2 the first Mott-lobe is the last to disappear. 
For all disorder strengths there is only one connected 
BG region, respectively one SF region. This is different 
from the system with tunneling disorder, but analogous 
to the on-site disordered case. 

A new feature occurs below the first Mott-lobe, which 
is shown as a blow-up in Figure 3d): In this region the 



FIG. 3: Phase diagram for box distributed disordered 
two-particle interactions Ui with the disorder strength 
A/U — 0.35, 0.6, 0.95 and a blow-up for A = 0.6 showing 
the tricritical point at p/U — JZ/U — 0. The black line is 
the perturbative result for MI-SF transition in the ordered 
system given by (2). 


BG-SF transition widely follows the MI-SF transition of 
the ordered case. Between this transition line and the 
MI-BG transition line the BG phase forms a narrowing 
strip. Both transition lines approach each other tightly 
for small tunneling rates and form a tricritical point 
in the limit fi/U —)► 0 and JZ/U 0, which does not 
contradict the fact that a direct transition from MI to 
SF is impossible in the disordered case [9, 38]. 

This phenomenon can be understood recalling 
equation (2) and studying the behavior of this equation 
under a variation of U : From its derivation, equation (2) 
describes the transition line, at which the SF order 
parameter ^ becomes non-zero in the ordered case. 
In the disordered case this takes place at the BG-SF 
transition. Therefore, the BG-SF transition line in the 
disordered case follows equation (2) in regions, where 
it is stable against variation of U. This variation of 
/jL± (J, Z, £/, n), which results from perturbation theory 
for small J [3], under a change of U is given by the 
derivative The expansion of the derivative for 

small J is given by 


d/i± 

dU 


x 2 n (n + 1) + n upper branch, 

— x 2 n (n + 1) + n — 1 lower branch, 


(7) 


in powers of x = Notice that the linear term cancels 
and in general is different from zero. Only in case of the 
lower branch of the first Mott-lobe (n = 1) it vanishes 
for zero tunneling rate. For increasing tunneling rates 
it grows less then linearly, since x is smaller than one. 
This means that the lower branch of the first Mott-lobe 
is fairly stable against variation of U. For all other Mott- 
lobes n > 1 the absolute value of the derivative is posi¬ 
tive for small tunneling rates J. This feature of the lower 
branch of the first Mott-lobe is unique and does not occur 
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for other branches of the disordered inter-particle inter¬ 
action case. Therefore, the BG-SF transition of the in¬ 
ter particle interaction disordered system below the first 
Mott-lob widely follows the transition line of the ordered 
system, which is given by equation ( 2 ), leading to the 
tricritical point at the origin of the phase diagram. 

IV. BICHROMATIC POTENTIAL 

Experimentally disorder can be introduced either by a 
diffuser [31] or by a bichromatic potential [4, 31, 32]. The 
diffuser modifies the intensity of the laser, which leads to 
inhomogeneities in the resulting optical lattice. For a de¬ 
tailed comparison with theoretical predictions a thorough 
characterization of the diffuser is necessary. Especially 
the width of the disorder distribution is a crucial system 
parameter, which is fixed by the diffuser and cannot be 
tuned freely. Alternatively a quasi-periodic potential is 
formed by a main optical lattices with a high intensity, 
which is superposed by a second weaker one with slightly 
different wave length [32]. The resulting lattice is not 
periodic but quasi-periodic and displays local inhomo¬ 
geneities. Such a quasi-periodic potential is the basis for 
our calculation, from which we will extract all BH pa¬ 
rameters and finally discuss the resulting phase diagram 
in dependence of the laser intensities of both lattices for 
integer filling. 

The quasi-periodic potential in two dimensions is given 
by 

V ( x , z) = Vl (cos 2 (kix) + cos 2 (kiz)) 

-hV 2 (cos 2 (k 2 x) + cos 2 (k 2 z)) , ( 8 ) 

where ki = y 1 (z = 1 , 2 ), the lattice constant a = ^ and 
the intensities V = s^Eri are given in units of the recoil 
energy Eri = • The wavelengths are chosen to be 

Ai = 830 nm and Ai = 1076 nm with reference to the 
experimental setup of [4]. In experiments 87 Rb, which 
has a mass of m m 1.44310 -25 kg, is widely used. The 
amplitude of the main lattice si determines the depth of 
the lattice. The amplitude of the second lattice s 2 <C si 
is by far smaller than the first one and increases the in¬ 
fluence of the disorder strength. 



FIG. 4: A quasi-periodic potential in one dimension according 
to V ( x) — Vi cos 2 (kix) + V 2 cos 2 (k 2 x) with si = 16, s 2 — 2.5 
and V 2 /Vi w0.09. 


A. Wannier-Functions 

In order to map the quasi-periodic potential ( 8 ) to the 
BH Hamiltonian the Wannier functions for each site are 
computed individually. We first focus on a regular pe¬ 
riodic (V 2 = 0) lattice. According to the Bloch theorem 
the Bloch functions 

r n (f) = u q -(r) = y>f 5 e- i( ^ (9) 

G 

solve the stationary Schrodinger equation. For ev¬ 
ery wave-vector k there exists a unique decomposition 
k = q — G, where q lives in the first Brillouin zone (l.BZ). 
The Bloch coefficients u^(r) are periodic functions with 
the same periodicity as the lattice. The Bloch functions, 
which spread over the hole lattice, form an orthonor¬ 
mal basis. Thus, the Wannier functions localized at site 
i = (i x ,i z ) can be construed as 

W n^ = J~h E xi=la. (10) 

V 1 qel.BZ 

They are real functions, which fulfill f y dVW^ 2 (r) = 1. 
Moreover, they are symmetric due to the underlying lat¬ 
tice symmetry. In Figure 5 a) a Wannier function for a 
regular ID lattice is shown exemplaryly. This wave func¬ 
tion shows a dominant occupation probability at one sin¬ 
gle site and small probability at the neighboring sites. 

In the non-symmetric case (0 7 ^ V 2 <C V±) the functions 

= ai) 

G 

still form an orthonormal basis, but their coefficients 
u$ (r) are no longer periodic. But still, localized func¬ 
tions according to generalized Wannier functions can be 
constructed according to equation ( 10 ). As one can see 
in Figure 5 on the right, these functions are still local¬ 
ized at a specific lattice site, but they are asymmetric, 
reflecting the asymmetry of the underlying lattice. 

With this generalization of the Wannier function for 




FIG. 5: Wannier functions for an ordered symmetric lattice 
on the left (s\ = 2, s 2 = 0) and a disordered asymmetric 
lattice on the right (s 1 = 2, s 2 = 1 with V 2 /Vi ~ 0.3). 

asymmetric lattice systems we follow the common deriva¬ 
tion introducing the BH model in literature [39, 40] and 
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thus fundamentally introduce disorder via quasi-periodic 
potentials. Other approaches avoid this step and describe 
this effect on the BH parameter effectively [15, 16, 41? 
, 42]. As we were able to show it is not sufficient to use 
the symmetric Wannier functions [44, 45] as an approx¬ 
imation in the disordered case; this is only accurate for 
52 < 0.1. Thus, in the following we will use these gen¬ 
eralized Wannier functions in order to determine the BH 
parameters. 


which especially the tunneling rate is negative. Finally, 
the chemical potential /i as a Lagrange multiplier for the 
condition of a fixed particle number N = JT(rq) = M of 
on average one particle per site is determined recursively 
within LMF theory. 

2. Distributions 


B. Bose-Hubbard Parameter 


The result of this procedure is a set of individual BH pa- 


1. Determination 


Bosons in any potential V (r) in the quasi-ideal regime 
are described by the quantum field theory Hamiltonian 

H = J df # (r) 

+ i J dr dr" ’3> t (r) ^ (f) U$ {?) '5 (r) , (12) 

where U (r, r*) describes the two particle interaction [39]. 
The field operator in tight-binding approximation 


fi 2 V 2 
2 m 


+ V (r) U (f) 


$(r ) = J2w i (r)a i (13) 

i 


can be composed by the Wannier functions 
Wi (r) = Wq (f) of the lowest band (n = 0) and 

the creation operator a i creating particle at site 
i = (i x , iz)- In the tight-binding approximation the 
inter-particle interaction reduces to a point-like in¬ 
teraction U (r,r f ) = UoS(r — r f ) [40]. The effective 
inter-particle interaction in 2d is given by 


h 2 a s 18m7TUJ z 

m V h 


5.56 10 -11 h, 


(14) 


where a s = 5.2 nm is the scattering length, 

m = 1.44310 -25 kg is the mass of the 87 Rb atoms and 
uj z = 67 r kHz is the frequency of the vertical confinement 
[46-48]. The BH parameters may be calculated using the 
ground state Wannier function Wi (r) = Wq (r) and the 
actual potential V (f): 




Ui 


Jin 


J d^,(rl(-^ + V(rl) 

U 0 j d fWi 4 (r) 

J dr Wj (r) (— + V (r)^j 


Wi(r) 


Wi{r) 


(15) 


In contrast to the symmetric potential, all these integrals 
are not necessarily positive in the asymmetric case. Here 
the generalized Wannier functions (11) as well as the po¬ 
tential are asymmetric and rarely configurations occur, in 





FIG. 6: BH Parameter for a quasi-periodic potential accord¬ 
ing to equation (8) with si = 10 and S 2 = 0.5. The tunneling 
rate at site i — ( i x ,i z ) gives the value for the tunneling rate 
to the neighbor j = (i x + Note that this contains all 

information, since the tunneling rate Jij (15) is symmetric 
under a change of the indices and the Potential V (x, z) (8) 
under a change of the coordinates x and z. 


rameters a = e, J, U for each fixed parameter set (si, 52 ). 
Exemplary the BH parameters for si = 10 and 52 = 0.5 
are shown in Figure 6. In the figures each pixel rep¬ 
resents the value of e*, J^-, Ui at a specific site. The 
BH parameter follow the modulation of the lattice po¬ 
tential. Therefore, we now deal with distributions P (a), 
which depend on the parameter set ( 51 , 52 ) chosen for 
the amplitudes of the lasers. We will especially focus 
on their mean value a = f da aP (a) and the variance 
a & = f d aa 2 P(a). Since we want to compare results 
to the box distributed case from section III, where the 

variance is given by cr 2 = -jf, we define the width of 
the distribution A a according to this equation. With 
the help of both benchmarks we are able to compare the 
distributions with the scenarios of disorder in only one 
BH parameter, as introduced in section III. The mean 
value of the distribution P (a) here matches the site in¬ 
dependent BH parameters a = e, J, U from secti on III, 
while the width of the distribution = ^12(7^ corre¬ 
sponds to the disorder strength given as free parameter 
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FIG. 7: The mean value a (first row) and the width A a = \Zl2a%, (second row) with the variance of the distributions of 
the BH-Parameter a = e, J, U in units of the recoil energy Er±. Notice that in the figures of e and U the s i- and S 2 -axis are 
switched in comparison to the other figures. 



A c e /U = l A c j/U = « 8.5810“ 2 A c v /U = 2 


FIG. 8: The width A a = y/12 oj of the distributions of the BH-Parameter a = e, J, U in units U. The white line shows the 
the critical width A}, where the last Mott-lobe disappears, which is stated below the plot. 


in section III. 

The resulting mean value and width of the distribution 
are shown in Figure 7. The amplitude of the main lat¬ 
tice si is one order of magnitude larger than the one of 
the second lattice S 2 - In a shallow lattice (si small) the 
mean value of the on-site energy e and the inter-particle 
interaction U are small and grow with increasing depth 
of the lattice (si large). The mean value of the tunnel¬ 
ing rate J reaches its maximal value in a shallow lattice 
and decreases in a deep lattice. All mean values are in¬ 
dependent of the strength S 2 of the second lattice. The 
width of the distribution of the on-site energy A e is in¬ 
dependent of the amplitude s\ of the main lattice, but 
increases with the amplitude S 2 of the second one. As ex¬ 
pected, the amplitude S 2 of the second lattice indeed in¬ 
creases the disorder strength in the system. The width of 
the distributions of the tunneling rate A j and the inter¬ 
particle interaction Ajy show similar behavior, depending 
on both parameters s± and S 2 , however, their maximal 
values differ substantially. Both show increasing widths 
for increasing S 2 and adopt the maximal values for a shal¬ 
low (si small) and strongly disordered (s 2 large) lattice. 
The critical disorder strength at which the last Mott-lobe 
disappears, is A C e /U = 1 for pure on-site energy disorder 
and A j/U = 3 ~ 2 c ^ 2 « 0.0858 for pure tunneling disor¬ 
der. Above these vales only the BG and SF phase re¬ 


main. A comparison of the width of the distribution in 
units of V (see Figure 8) with the results on disorder in 
only one BH parameter in section III shows that both, 
the width of the distribution of the on-site energy, as well 
as the tunneling rate, reaches the region where all three 
phases occur in the phase diagram. Even though the oc¬ 
curring width of the distribution of the tunneling rate 
Aj is small, it reaches the parameter range, where all 
three phases compete in the phase diagram. In contrast, 
the width of the inter-particle interaction A*y is indeed 
small in comparison to the range in which all three phases 
occur in the phase diagram and thus may be neglected. 


V. PHASE DIAGRAMS OF THE 
BICHROMATIC POTENTIAL 

A. Intensity phase diagrams 


With the help of the determined BH parameters we have 
determined the phase diagram for fixed average number 
N = YU&i) = Af of one particle per site in dependence 
of the laser intensities s\ and si by adjusting the chemi¬ 
cal potential /i. The resulting phase diagram is shown in 
Figure 9 a). Here the BH parameters are correlated ac- 
















cording to the lattice potential (8), exemplarily displayed 
in Figure 6. According to the LMF cluster analysis de¬ 
scribed in section II, the MI phase is characterized by 
the absence of any SF site, which means that every site 
in the MI region has an integer particle number. In the 
SF region sites with non-integer particle number perco¬ 
late. In between, in the BG region, the system consists 
of both, sites with integer, and others with non-integer 
particle number, which do not percolate. 

Next we compare the obtained phase diagram for the 
bichromatic quasi-periodic potential with the one ob¬ 
tained for a BH model with uncorrelated disorder ac¬ 
cording to identical distributions of the BH parameters 
P (ol) with a = e, J, U. These distributions depend on 
both of the laser intensities si and 52. We start with 
one parameter set for e,J and U given by fixed 81 and 
s 2, which we have determined and discussed in sections 
IV B 1 and IV B 2. We produce 200 different samples, by 
randomly choosing new site indices. In other words, we 
study 200 samples according to the same distribution by 
switching lattice sites thereby erasing local correlations 
in the parameter set. After performing the LMF cluster 
analysis we determine the BG-SF boundary with finite 
size scaling. The resulting phase diagram in Figure 9 b) 
does not differ much from Figure 9 a), only the BG-SF 
transition line for uncorrelated disorder is slightly dis¬ 
torted in comparison with the quasi-periodic case. 

In the resulting phase diagram, shown in Figure 9 b), all 
three phases occur in dependence of the lattice parame¬ 
ters si and <§2• Along the si-axis (S2 = 0) the direct SF- 
MI transition of the ordered system occurs. For values 
below this point in a shallow lattice the SF phase cov¬ 
ers the whole parameter region independent of S2. This 
corresponds to the fact that in this region the tunneling 
rate is largest, as shown in Figure 7. Above this point 
the MI occurs, which is completely surrounded by the 
BG for intermediate 82, which in turn is enclosed by the 
SF phase for even larger amplitudes of the second lattice 
82. Notice that the potential (8) reduces to the ordered 
case for 82 = 0 as well as for 81 = 0. Therefore, along the 
82-axis the system also undergoes a direct MI-SF transi¬ 
tion. In the region where 81 82 the second lattice is 

dominant and a similar structure occurs. 


B. BH parameter phase diagrams 


Let us now have a look at the phase diagram in depen¬ 
dence of the BH parameters. In the case of only one 
disordered BH parameter two possible representations of 
the phase diagram are common: For the first one [6, 26] 
the particle number is fixed to one particle per site, which 
fixes the chemical potential fi. Then, the phase transi¬ 
tions are shown in dependence of the inter-particle inter¬ 
action U/J and the disorder strength A/2 J. The second 
representation [21-23, 26] shows the Mott-lobes in de¬ 
pendence of the tunneling rate JZ/U and the chemical 



FIG. 9: Left: Phase diagram for the bichromatic quasi- 
periodic potential (8) for particle density (h-) = 1 in depen¬ 
dence of the laser intensities s 1 and S 2 for a system with 32 x 32 
sites. Right: The same phase diagram as on the left side but 
for a BH model with uncorrelated on-site, hopping strengths 
and interaction energy disorder with the distribution param¬ 
eters determined for the bichromatic quasi-periodic potential 
in section IV B 2. 


potential n/U for fixed disorder strength A/U. In this 
section we will discuss our results in both representations, 
keeping in mind that all BH parameters as well as all dis¬ 
order strengths are functions of 81 and 82, and thus are 
not independent of each other. 

The data from the (81,82)-phase diagram shown in Fig¬ 
ure 9 can be translated into a diagram similar to the 
first representation: According to Figure 7 the mean 
parameters e, J, U as well as the widths A e , A j, A u 
are functions of the two amplitudes 81 and 82. Since 
A u is two orders of magnitude smaller than the other 
widths, the inter-particle interaction U can be treated 
as a sharp value to a good approximation. As a re¬ 
sult the phase diagram can be visualized as the surface 
(U/J, A e /2J, Aj/2J) in three dimensions. This is shown 
in Figure 10, where each phase is colored differently. No¬ 
tice that with a quasi-periodic potential (8), which de¬ 
pends on the two amplitudes 81 and 82, only this sur¬ 
face in the BH parameter space can be reached, since all 
BH parameters are functions of 81 and 82 and dependent 
on each other. As a consequence disorder, e.g, where 
only one parameter is disordered while the others are 
fixed, cannot be reached in the phase diagram. Either it 
is an ordered (A e = A j = 0) or a completely disordered 
(A e 7^ 0, Aj ^ 0) system. This has two important im¬ 
plications: With a quasi-periodic potential neither the 
whole parameter space nor a pure on-site disorder can 
be realized. 

The phase diagram in Figure 10 shows all three phases: 
The BG phase (blue) separates the MI (dark gray) 
phase at strong inter-particle interaction U/J from the 
SF regime (light gray) at weak inter-particle interaction 
U/J. In Figure 10 b) the same phase diagram is show as 
a projection on the (U/J, A e /2J)-plane. 

The phase boundaries differ substantially from those of 
the BH model with uncorrelated disorder exclusively in 
the on-site energies [6, 26] due to the additional presence 
of disorder in the hopping strengths. More drastic is the 
difference between Figure 10 and the V 2 /J-U/J phase 
diagram predicted for a one-dimensional BH-model with 
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FIG. 10: Phase diagram for a quasi-periodic potential (8) in 
dependence of the BH Parameter. 


bichromatic quasi-periodicity only in the on-site poten¬ 
tial [15]. Since according to Figure 7 d) A e is propor¬ 
tional to V 2 = S 2 Er 2 , the phase diagram in Figure 10 
is directly comparable to Figure 1 (left) of [15] and Fig¬ 
ure 3 (bottom) of [16], both of which show that a direct 
MI-SF transition occurs. The latter is absent in Figure 
10, where an intervening BG phase occurs between the 
MI and SF phase. This might be explained by the fact 
that already a small amount of disorder in the hopping 
strengths strongly enlarges the BG regions in the phase 
diagram, as Figure 2 demonstrates. 

For the second representation we fix the weaker ampli¬ 
tude 52, which introduces disorder to the system, and 
study the system in dependence of si and fi. Since the 
tunneling rate J is a unique function of s\ and indepen¬ 
dent of 52, as shown in Figure 11 b), the 51-axis can 
easily be converted into a J-axis. In theoretical works 
disorder is usually introduced by bounded distributions 
with zero mean values. In the quasi-random case the 
mean value of the distributions P ( a ) of the BH- param¬ 
eters a = e, J, U are non-zero, as shown in Figure 7. In 
order to take this into account, we use fi — e instead of 
simply fi. Thus, from the data in the (51, fi)- plane, we 
can extract a phase diagram in dependence of JZ/U and 
(/i — e) /[/, as shown in Figure 12. Notice that J, [/, e as 
well as Aj, A*/, A e are all functions of s\ and 52 and not 
independent from each other, as they are asumed to be 
in most simulations of disordered systems. 

In Figure 11 the behavior of the system parameters for 
different values of 5 2 = 0.0354, 0.0758, 0.1162 is shown as 
a function of JZ/U. Figure b) shows the tunneling rate in 
dependence of s\. While the tunneling rate JZ/U is inde¬ 
pendent of 52, the on-site energy e/U varies for different 
values of 52. In a shallow lattice (51 small) the tunneling 
rate JZ/U is large, while for increasing s\ it approaches 
zero and finally in a deep lattice (51 large) the tunnel¬ 
ing rate JZ/U becomes infinitesimally small. This means 
that in the phase diagram the fi/U -axis at JZ/U = 0 may 
be approached with arbitrary accuracy, but can never be 
reached. In a deep lattice ( JZ/U small, s\ large) the 
on-site energy e/U is small and increases with growing 
JZ/U. The disorder strengths A e /U, Aj/U , A u/U in¬ 
crease with JZ/U and the amplitude 52. 

The phase diagram for different values of 52 as a function 
of JZ/U and (fi — e) /U is shown in Figure 12. Notice 



FIG. 11: The on-site energy e/U and the widths 

A 6 /U : A j/Uj Au/U of the BH parameter for different values 

of s 2 = 0.0354 (-) , 0.0758 (-•) , 0.1162 ( - ) as a function of 

JZ/U. In the middle on the top the tunneling rate JZ/U is 
shown as a function of si. The tunneling rate is independent 
of S 2 and a unique function of s±. The on-site energy e/U as 
well as all widths increase with the tunneling rate JZ/U and 
the amplitude S 2 . 


that the BH parameters and their widths vary along the 
JZ /[/-axis corresponding to Figure 11. For all values of 
52 we find a regular structure of Mott-lobes, surrounded 
by individual BG regions. The number of Mott-lobes as 
well as the number of BG regions decreases with increas¬ 
ing disorder amplitude 52, which is equal to the increase 
of the disorder strengths A e /C7, Aj/C7, A u/U. The reg¬ 
ular pattern of Mott-lobes and BG regions repeats it¬ 
self in intervals of length one along the (fi — e ) /[7-axis. 
The lower and the upper extent of the Mott-lobes have 
the same distance to the next integer number for fixed 
52. Thus, the Mott-lobes have the same width along the 
(fi — e) /[/-axis, while their extension in JZ/U -direction 
shrinks with their number n. Except of the first BG re¬ 
gion all the others are separated from each other by SF re¬ 
gions, reaching down to very small tunneling rates JZ/U. 
This is a unique feature of disorder only in the tunneling 
rates, which was discussed in section IIIB. The fact that 
we see this special phenomenon here in the phase dia¬ 
gram of a quasi-random potential, once more promotes 
our finding from section IV B 2 that the influence of dis¬ 
order in the tunneling rate cannot be neglected. 


VI. CONCLUSION 

The main result of this paper is that disorder in the 
tunneling rate must be taken into account in setups using 
2d bichromatic quasi-periodic potentials. It influences 
the phase diagram just as much as the on-site disorder. 
This claim is supported by several findings in this paper: 
Firstly, in section III we showed that the disorder 
strength, where all three phases compete in the phase 
diagram, is one order of magnitude smaller for tunneling 
disorder than for on-site interaction. We found that each 
scenario, in which only one BH parameter is disordered, 
yields different features in the phase diagram. Especially, 
the characteristics of tunneling disorder are a finite num- 
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FIG. 12: Phase diagram for different 

S 2 = 0.0354, 0.0758, 0.1162 as a function of the mean 
tunneling rate J Z/U and the chemical potential (fi — e) /U. 


ber of Mott-lobes and the existence of SF regions even 
for JZ/U = 0. 

Secondly, we discussed bichromatic quasi-periodic poten¬ 
tials in section IV and showed (see Figure 7) that the 
width of the distribution of the tunneling rate as well as 
that of the on-site energy reach the physical interesting 
region, where all three phases compete. This is true even 
though the width of the distribution of the tunneling rate 
is one order of magnitude smaller than that of the on-site 
energy. The influence of disordered inter-particle inter¬ 
action is negligible, since its width remains four orders 
of magnitude below the critical disorder strength. This 
is in agreement with the results for distributions of BH 
parameters produced by a random diffuser pattern over¬ 
lapping the main lattice [31, 33]. Correspondingly to our 
work, in these papers the authors also showed, that the 
width of the tunneling rate and the inter-particle inter¬ 
action are several orders of magnitude smaller than that 
of the on-site energy. Moreover, it is important to keep 
in mind that in a bichromatic quasi-periodic potential it 
is not possible to study exclusively on-site disorder. For 
growing intensity of the second laser S2, the widths of 
the disorder distributions of both, the on-site, and the 
tunneling rate, increase simultaneously. 

Thirdly, the influence of tunneling disorder is obvious in 
the phase diagram of the quasi-periodic potential in de¬ 
pendence of the BH parameters. The transition lines, 
shown in Figure 10, deviate from pure box distributed 
on-site disorder, discussed for example in [6, 26]. In the 
quasi-random case the SF region is smaller, while the 
BG and the MI regions cover a larger region. In the 
fi-J- phase diagram of Figure 12 we find individual BG 
regions, which are separated by SF regions, which is a 


unique feature exclusively occurring in systems with dis¬ 
ordered tunneling rates (see Figure 2). 

While the field of box on-site disorder was studied widely 
[6-26, 35], the works on tunneling disorder are rare and 
mainly deal with bimodal distributions [10, 36, 37]. A 
phase diagram far above the critical disorder strength 
A j/U ~ 8.58 10 -2 for equally distributed tunneling dis¬ 
order is shown in [25]. 

The one-dimensional BH model with a bichromatic on¬ 
site potential with incommensurable wave lengths was 
studied in [15, 16]. The motivation there was, as also in 
this paper, to qualitatively understand the phase diagram 
of bosons in a bichromatic quasi-periodic potential, but it 
was argued that the variations in the hopping strengths 
as well as in the in the interaction energies were only mi¬ 
nor and could be neglected. Although, we have shown 
here by explicit calculation that the disorder strength 
in the tunneling strengths is indeed one order of mag¬ 
nitude smaller than the on-site disorder, it nevertheless 
has a strong effect on the 2d phase diagram. The most 
striking difference is that the phase diagrams of [15, 16] 
show a direct MI-SF transition (which for quasi-period 
disorder is not in contradiction with general predictions 
for uncorrelated disorder [9]), whereas we find an inter¬ 
vening BG phase between the MI and the SF phase for 
all values of the laser intensities s i and 52- A reason 
for this, in addition to potentially qualitative differences 
between one and two dimensions, could be that our re¬ 
sults for uncorrelated disorder in section III show that a 
modest amount of disorder in the hopping strength al¬ 
ready generates relatively large BG regions in the phase 
diagram (see Figure 2). Consequently, disorder in the 
hopping strengths cannot be neglected studying bosons 
in a bichromatic quasi-periodic potential. 

Conversely, one should be aware that experimental re¬ 
alizations of the disordered potential by a bichromatic 
quasi-period potential, as in [4, 32], produce a phase dia¬ 
gram that is qualitatively very different from the predic¬ 
tions of the disordered BH model with exclusively on-site 
disorder. 
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